program example_FEM
	use MysparseOperator
	use myFEM
	use common_functions
	implicit none
	double precision :: domain(4), O(2), width, height
	integer :: m, n

	double precision, allocatable :: u(:)
	double precision :: errors(4)


	integer :: dofs, ii, iter, tt


	! the domain data
	O = (/ 0.d0, 0.d0 /)
	width = 1.d0
	height = 1.d0

	write(*,'("input the number of refinements "$)')
	read(*,*) iter
	write(*,*) 'input the initial mesh data:'
	write(*,'("      m0 = "$)')
	read(*,*) m
	write(*,'("      n0 = "$)')
	read(*,*) n
! 	iter = 2
! 	m = 4
! 	n = 4


	do ii = 1, iter, 1
		dofs = (m+1)*(n+1)
		allocate(u(dofs))
		call mf_linear_cartesian_DB_dsolver(O,width,height,m,n,a_fun_d,c_fun_d,f_fun_d,g_fun_d,u)
		errors = 1.d0
		call mf_linear_cartesian_derror(O,width,height,m,n,u_fun_d,u,errors)
		write(*,'(I4,4E12.3)') m, errors
		m = m*2
		n = n*2

		if ( allocated(u) )  deallocate(u)
	end do


contains


	subroutine u_fun_d(x,y,v,pxy)
		double precision, dimension(:), intent(in) :: x, y
		double precision, dimension(:), intent(inout) :: v
		character, intent(in), optional :: pxy
		if ( .not.present(pxy) ) then
			v = x*(1.d0-x)*y*(1.d0-y)
		end if
		if ( present(pxy) .and. pxy.eq.'x' ) then
			v = y*(1.d0-y)*(1.d0-2.d0*x)
		end if
		if ( present(pxy) .and. pxy.eq.'y' ) then
			v = x*(1.d0-x)*(1.d0-2.d0*y)
		end if
	end subroutine u_fun_d

	subroutine a_fun_d(x,y,v)
		double precision, dimension(:), intent(in) :: x, y
		double precision, dimension(:), intent(inout) :: v
		v = 1.d0
	end subroutine a_fun_d

	subroutine c_fun_d(x,y,v)
		double precision, dimension(:), intent(in) :: x, y
		double precision, dimension(:), intent(inout) :: v
		v = 1.d0
	end subroutine c_fun_d
	subroutine f_fun_d(x,y,v)
		double precision, dimension(:), intent(in) :: x, y
		double precision, dimension(:), intent(inout) :: v
		v = 2.d0*x*(1-x) + 2.d0*y*(1-y) + x*(1-x)*y*(1-y)
	end subroutine f_fun_d
	subroutine g_fun_d(x,y,v)
		double precision, dimension(:), intent(in) :: x, y
		double precision, dimension(:), intent(inout) :: v
		v = 0.d0
	end subroutine g_fun_d




end program example_FEM
